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Abstract 



Given i.i.d. observations of a random vector X G R*", we study the problem of estimating both its covariance 
matrix E*, and its inverse covariance or concentration matrix B* = (E*)~^. We estimate O* by minimizing an 
■ ^1 -penalized log-determinant Bregman divergence; in the multivariate Gaussian case, this approach corresponds to 

^i-penalized maximum likelihood, and the structure of O* is specified by the graph of an associated Gaussian Markov 
random field. We analyze the performance of this estimator under high-dimensional scaling, in which the number of 
nodes in the graph p, the number of edges s and the maximum node degree d, are allowed to grow as a function of the 
sample size n. In addition to the parameters (p, s, d), our analysis identifies other key quantities that control rates: 
C/3 . (a) the i?oo -operator norm of the true covariance matrix E*; and (b) the £oo operator norm of the sub-matrix Fgs, 

where S indexes the graph edges, and F* = (0*)~^ (g) and (c) a mutual incoherence or irrepresentability 

measure on the matrix F* and (d) the rate of decay l//(n, 5) on the probabilities {|E^ — E*j| > 5}, where E" 
^ ' is the sample covariance based on n samples. Our first result establishes consistency of our estimate O in the ele- 

00 , mentwise maximum-norm. This in turn allows us to derive convergence rates in Frobenius and spectral norms, with 

improvements upon existing results for graphs with maximum node degrees d = o{^/s). In our second result, we 
show that with probability converging to one, the estimate O correctly specifies the zero pattern of the concentration 
matrix O*. We illustrate our theoretical results via simulations for various graphs and problem parameters, showing 
good correspondences between the theoretical predictions and behavior in simulations. 



OO 

9^ ; 1 Introduction 

\ The area of high-dimensional statistics deals with estimation in the "large p, small n" setting, where p and n corre- 

■ spond, respectively, to the dimensionality of the data and the sample size. Such high-dimensional problems arise in a 

I variety of applications, among them remote sensing, computational biology and natural language processing, where 

the model dimension may be comparable or substantially larger than the sample size. It is well-known that such high- 
dimensional scaling can lead to dramatic breakdowns in many classical procedures. In the absence of additional model 
assumptions, it is frequently impossible to obtain consistent procedures when p ^ n. Accordingly, an active line of 
statistical research is based on imposing various restrictions on the model — for instance, sparsity, manifold structure, 
or graphical model structure — and then studying the scaling behavior of different estimators as a function of sample 
size n, ambient dimension p and additional parameters related to these structural assumptions. 

In this paper, we study the following problem: given n i.i.d. observations {X'-'^^^^i of a zero mean random vector 
X e RP, estimate both its covariance matrix E*, and its inverse covariance or concentration matrix 6* := (S*) ^. 
Perhaps the most natural candidate for estimating E* is the empirical sample covariance matrix, but this is known to 
behave poorly in high-dimensional settings. For instance, when p/n ^ c > Q, and the samples are drawn i.i.d. from 
a multivariate Gaussian distribution, neither the eigenvalues nor the eigenvectors of the sample covariance matrix 
are consistent estimators of the population versions [14, 15]. Accordingly, many regularized estimators have been 
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proposed to estimate the covariance or concentration matrix under various model assumptions. One natural model 
assumption is that reflected in shrinkage estimators, such as in the work of Ledoit and Wolf [ 1 6], who proposed 
to shrink the sample covariance to the identity matrix. An alternative model assumption, relevant in particular for 
time series data, is that the covariance or concentration matrix is banded, meaning that the entries decay based on 
their distance from the diagonal. Furrer and Bengtsson [II] proposed to shrink the covariance entries based on this 
distance from the diagonal. Wu and Pourahmadi [24] and Huang et al. [13] estimate these banded concentration 
matrices by using thresholding and -penalties respectively, as applied to a Cholesky factor of the inverse covariance 

matrix. Bickel and Levina [!'] prove the consistency of these banded estimators so long as ^^2SJ^ y q and the model 

covariance matrix is banded as well, but as they note, these estimators depend on the presented order of the variables. 

A related class of models are based on positing some kind of sparsity, either in the covariance matrix, or in the 
inverse covariance. Bickel and Levina [ i ] study thresholding estimators of covariance matrices, assuming that each 
row satisfies an £q-ball sparsity assumption. In independent work. El Karoui [9] also studied thresholding estimators 
of the covariance, but based on an alternative notion of sparsity, one which captures the number of closed paths of any 
length in the associated graph. Other work has studied models in which the inverse covariance or concentration matrix 
has a sparse structure. As will be clarified in the next section, when the random vector is multivariate Gaussian, the 
set of non-zero entries in the concentration matrix correspond to the set of edges in an associated Gaussian Markov 
random field (GMRF). In this setting, imposing sparsity on the concentration matrix can be interpreted as requiring 
that the graph underlying the GMRF have relatively few edges. A line of recent papers [8, 10, 25] have proposed an 
estimator that minimizes the Gaussian negative log-Ukelihood regularized by the £i norm of the entries (or the off- 
diagonal entries) of the concentration matrix. The resulting optimization problem is a log-determinant program, which 
can be solved in polynomial time with interior point methods [ j], or by faster co-ordinate descent algorithms [8, n ,]. In 
recent work, Rothman et al. [21] have analyzed some aspects of high-dimensional behavior of this estimator; assuming 
that the minimum and maximum eigenvalues of S* are bounded, they show that consistent estimates can be achieved 

in Frobenius and operator norm, in particular at the rate 0{^J^^^^). 

The focus of this paper is the problem of estimating the concentration matrix 0* under sparsity conditions. We 
do not impose specific distributional assumptions on X itself, but rather analyze the estimator in terms of the tail 
behavior of the maximum deviation max^ j 1 1]"^ — S*j | of the sample and population covariance matrices. To estimate 
Q*, we consider minimization of an £i-penalized log-determinant Bregman divergence, which is equivalent to the 
usual £i-penaUzed maximum likelihood when X is multivariate Gaussian. We analyze the behavior of this estimator 
under high-dimensional scaling, in which the number of nodes p in the graph, and the maximum node degree d are all 
allowed to grow as a function of the sample size n. 

In addition to the triple {n,p, d), we also explicitly keep track of certain other measures of model complexity, 
that could potentially scale as well. The first of these measures is the £oo-operator norm of the covariance matrix E*, 
which we denote by /\y;. := |||oo- The next quantity involves the Hessian of the log-determinant objective function, 
r* (O*)^^ (g) (9*)^^. When the distribution of X is multivariate Gaussian, this Hessian has the more explicit 
representation F^- ^.^ = cov{XjXk, X^Xm}, showing that it measures the covariances of the random variables 
associated with each edge of the graph. For this reason, the matrix F* can be viewed as an edge-based counterpart to 
the usual node-based covariance matrix E*. Using S to index the variable pairs associated with non-zero entries 
in the inverse covariance. our analysis involves the quantity Kr* = lll(r55)~^|||oo- Finally, we also impose a mutual 
incoherence or irrepresentability condition on the Hessian F*; this condition is similar to assumptions imposed on E* 
in previous work [22, 26, 19, 23] on the Lasso. We provide some examples where the Lasso irrepresentability condition 
holds, but our corresponding condition on F* fails; however, we do not know currently whether one condition strictly 
dominates the other. 

Our first result establishes consistency of our estimator 8 in the elementwise maximum-norm, providing a rate 
that depends on the tail behavior of the entries in the random matrix E" — E*. For the special case of sub-Gaussian 
random vectors with concentration matrices having at most d non-zeros per row, a corollary of our analysis is consis- 
tency in spectral norm at rate |||0 — &* |||2 = {d'^ ^ogp)/n), with high probabihty, thereby strengthening previous 
results [21]. Under the milder restriction of each element of X having bounded 4m-th moment, the rate in spectral 
norm is substantially slower — namely, |||8 — 8* |||2 — 0{dp^^^^'^ / ^/n) — highlighting that the familiar logarithmic de- 
pendence on the model size p is linked to particular tail behavior of the distribution of X. Finally, we show that under 
the same scalings as above, with probability converging to one, the estimate 8 correctly specifies the zero pattern of 
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the concentration matrix Q*. 

The remainder of this paper is organized as follows. In Section 2, we set up the problem and give some background. 
Section 3 is devoted to statements of our main results, as well as discussion of their consequences. Section 4 provides 
an outline of the proofs, with the more technical details deferred to appendices. In Section 5, we report the results of 
some simulation studies that illustrate our theoretical predictions. 

Notation For the convenience of the reader, we summarize here notation to be used throughout the paper Given a 
vector M G M'' and parameter a G [1, oo], we use ||w|jQ to denote the usual norm. Given a matrix U G MP^p and 
parameters a,b G [1, oo], we use |||?7|||a.b to denote the induced matrix-operator norm max||j,||^^i see Horn 

and Johnson [12] for background. Three cases of particular importance in this paper are the spectral norm |||C/|||2, 
corresponding to the maximal singular value of U; the ioo/ ^oa-operator norm, given by 

p 

\\U\\oo .max ^|C/,fc|, (1) 

and the £i / £i-operator norm, given by |||[/|||i = ||| C^"'" ||| oo ■ Finally, we use j|C/j|oo to denote the element-wise maximum 
maxi j \Uij\; note that this is not a matrix norm, but rather a norm on the vectorized form of the matrix. For any 
matrix U G R^^^, we use vec(C/) or equivalently U £ W to denote its vectorized form, obtained by stacking 
up the rows of U. We use {{U, V)) := J2i j ^ij^ij ^'^ denote the trace inner product on the space of symmetric 

matrices. Note that this inner product induces the Fro^en;Mi noniz |||C/|||f := \J^^i~jUfj- Finally, for asymptotics, 

we use the following standard notation: we write J{n) = 0{g{n)) if f{n) < cg{n) for some constant c < oo, and 
/(n) = ^{g{n)) if f{n) > c'g{n) for some constant c' > 0. The notation f{n) x g{n) means that f{n) = 0{g{n)) 
and/(n) = ^{g{n)). 



2 Background and problem set-up 

Let X — (Xi, . . . , Xp) be a zero mean p-dimensional random vector The focus of this paper is the problem of 
estimating the covariance matrix S* := E[XX-^] and concentration matrix 0* :— E*^^ of the random vector X 
given n i.i.d. observations {X'^^^lJJ^j. In this section, we provide background, and set up this problem more pre- 
cisely. We begin with background on Gaussian graphical models, which provide one motivation for the estimation of 
concentration matrices. We then describe an estimator based based on minimizing an ii regularized log-determinant 
divergence; when the data are drawn from a Gaussian graphical model, this estimator corresponds to -regularized 
maximum likelihood. We then discuss the distributional assumptions that we make in this paper. 

2.1 Gaussian graphical models 

One motivation for this paper is the problem of Gaussian graphical model selection. A graphical model or a Markov 
random field is a family of probability distributions for which the conditional independence and factorization properties 
are captured by a graph. Let X — {Xi,X2, . . . , Xp) denote a zero-mean Gaussian random vector; its density can be 
parameterized by the inverse covariance or concentration matrix Q* = (S*)^^ G S'^, and can be written as 

We can relate this Gaussian distribution of the random vector X to a graphical model as follows. Suppose we are 
given an undirected graph G — (V, E) with vertex set V = {1, 2, . . . ,p} and edge' set E, so that each variable Xi is 
associated with a corresponding vertex i G V. The Gaussian Markov random field (GMRF) associated with the graph 

' As a remark on notation, we would like to contrast the notation for the edge-set E from the notation for an expectation of a random variable, 
E(.). 
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Zero pattern of inverse covariance 




(a) 
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(b) 



Figure 1. (a) Simple undirected graph. A Gauss Markov random field has a Gaussian variable Xi associated with each 
vertex i £ V. This graph has p — 5 vertices, maximum degree d = 3 and s = 6 edges, (b) Zero pattern of the inverse 
covariance 0* associated with the GMRF in (a). The set £(0*) corresponds to the off -diagonal non-zeros (white blocks); 
the diagonal is also non-zero (grey squares), but these entries do not correspond to edges. The black squares correspond to 
non-edges, or zeros in 0*. 



G over the random vector X is then the family of Gaussian distributions with concentration matrices Q* that respect 
the edge structure of the graph, in the sense that 8*^ = if ^ E. Figure 1 illustrates this correspondence between 
the graph structure (panel (a)), and the sparsity pattern of the concentration matrix O* (panel (b)). The problem of 
estimating the entries of the concentration matrix Q* corresponds to estimating the Gaussian graphical model instance, 
while the problem of estimating the off-diagonal zero-pattern of the concentration matrix — that is, the set 

Eie*) := {i,jeV\ i^j,e*,^0} (3) 

corresponds to the problem of Gaussian graphical model selection. 

With a slight abuse of notation, we define the sparsity index s :~ \E{Q*) \ as the total number of non-zero elements 
in off-diagonal positions of Q*; equivalently, this corresponds to twice the number of edges in the case of a Gaussian 
graphical model. We also define the maximum degree or row cardinality 



max 
:=i,...,7 



{]eV\ % ^ 0} 



(4) 



corresponding to the maximum number of non-zeros in any row of O*; this corresponds to the maximum degree in the 
graph of the underlying Gaussian graphical model. Note that we have included the diagonal entry Q*^ in the degree 
count, corresponding to a self-loop at each vertex. 

It is convenient throughout the paper to use graphical terminology, such as degrees and edges, even though the 
distributional assumptions that we impose, as described in Section 2.3, are milder and hence apply even to distributions 
that are not Gaussian MRFs. 



2.2 £i -penalized log-determinant divergence 

An important set in this paper is the cone 

Si := {AeRP^P I ^^0}, (5) 

formed by all symmetric positive semi-definite matrices in p dimensions. We assume that the covariance matrix S* 
and concentration matrix 6* of the random vector X are strictly positive definite, and so lie in the interior of this cone 
SI. 

The focus of this paper is a particular type of 7\jf-estimator for the concentration matrix 6*, based on minimizing 
a Bregman divergence between symmetric matrices. A function is of Bregman type if it is strictly convex, continu- 
ously differentiable and has bounded level sets [4, 7]. Any such function induces a Bregman divergence of the form 
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Dg{A\\B) = g{A)~g{B) - {\/g{B),A- B). From the strict convexity of .g, it follows that Dg(A|lB) > for all A 
and B, with equality if and only if A — B. 

As a candidate Bregman function, consider the log-determinant barrier function, defined for any matrix A e 

by 

giA) :^ /-l«gdet(^) ifA^O 
I +00 otherwise. 

As is standard in convex analysis, we view this function as taking values in the extended reals R* = R U {+00}. 
With this definition, the function g is strictly convex, and its domain is the set of strictly positive definite matrices. 
Moreover, it is continuously differentiable over its domain, with Vg{A) = —A~^; see Boyd and Vandenberghe [3] for 
further discussion. The Bregman divergence corresponding to this log-determinant Bregman function g is given by 

Dg{A\\B) := -logdctA + logdctB + {{B-\A-B)), (7) 

valid for any A,B£ that are strictly positive definite. This divergence suggests a natural way to estimate concen- 
tration matrices — namely, by minimizing the divergence Dg (8* \\Q) — or equivalently, by minimizing the function 

rnin{((e, E*)) -logdcte}, (8) 

where we have discarded terms independent of &, and used the fact that the inverse of the concentration matrix is 
the covariance matrix (i.e., (9*)~^ = E* = E[XX^]). Of course, the convex program (8) cannot be solved without 
knowledge of the true covariance matrix E*, but one can take the standard approach of replacing E* with an empirical 
version, with the possible addition of a regularization term. 

In this paper, we analyze a particular instantiation of this strategy. Given n samples, we define the sample covari- 
ance matrix 

:= iyxW(xWf. (9) 

To lighten notation, we occasionally drop the superscript n, and simply write E for the sample covariance. We also 
define the off-diagonal £1 regularizer 

lieiii.off := ^|e,y|, (10) 

where the sum ranges over all i,j = 1, . . . ,p with i ^ j- Given some regularization constant A„ > 0, we consider 
estimating 8* by solving the following £i- regularized log-determinant program: 

e ~ argnmi{((e, E"))-logdct(e) + A„||e||i,off}. (11) 

As shown in Appendix A, for any A„ > and sample covariance matrix E" with strictly positive diagonal, this convex 
optimization problem has a unique optimum, so there is no ambiguity in equation (11). When the data is actually drawn 
from a multivariate Gaussian distribution, then the problem (1 1) is simply ^1 -regularized maximum likelihood. 



2.3 Tail conditions 

In this section, we describe the tail conditions that underlie our analysis. Since the estimator (11) is based on using 
the sample covariance E" as a surrogate for the (unknown) covariance E*, any type of consistency requires bounds on 
the difference E" — E*. In particular, we define the following tail condition: 

Definition 1 (Tail conditions). The random vector X satisfies tail condition T(/, u*) if there exists a constant v^, G 
(0, 00] and a function / : N x (0, 00) ^ (0, oo) such that for any (z, j) & V x V: 

P[|E^;. -E^l > J] < l//(n,<5) forallJe (0,1/w,]. (12) 

We adopt the convention 1/0 := +00, so that the value f * = indicates the inequality holds for any S G (0, 00). 
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Two important examples of the tail function / are the following: 

(a) an exponential-type tail function, meaning that f{n,5) = exp(c7i 5"), for some scalar c > 0, and exponent 
a > 0; and 

(b) a polynomial-type tail function, meaning that /(n, 5) = cti™ (5^™, for some positive integer to G N and scalar 
c> 0. 

As might be expected, if X is multivariate Gaussian, then the deviations of sample covariance matrix have an 
exponential-type tail function with a = 2. A bit more generally, in the following subsections, we provide broader 
classes of distributions whose sample covariance entries satisfy exponential and a polynomial tail bounds (see Lem- 
mata 1 and 2 respectively). 

Given a larger number of samples n, we expect the tail probability bound 1 / f{n, S) to be smaller, or equivalently, 
for the tail function f{n, S) to larger. Accordingly, we require that / is monotonically increasing in n, so that for each 
fixed (5 > 0, we can define the inverse function 

nf{r;S) :— argmaxjri | f{n,5)<r^. (13) 

Similarly, we expect that / is monotonically increasing in S, so that for each fixed n, we can define the inverse in the 
second argument 

6f{r;n) := argmax{(5 | f{n,S)<r^. (14) 
For future reference, we note a simple consequence of the monotonicity of the tail function / — namely 

n>nf{S,r) for some (5 > =^ Sf{n,r)<S. (15) 

The inverse functions fif and Sf play an important role in describing the behavior of our estimator We provide 
concrete examples in the following two subsections. 

2.3.1 Sub-Gaussian distributions 

In this subsection, we study the case of i.i.d. observations of sub-Gaussian random variables. 

Definition 2. A zero-mean random variable Z is sub-Gaussian if there exists a constant a E (0, oo) such that 

E[exp(tZ)] < exp(cr2iV2) foralHeM. (16) 

By the Chernoff bound, this upper bound (16) on the moment-generating function implies a two-sided tail bound 
of the form 

¥[\Z\>z] < 2cxp(-^). (17) 

Naturally, any zero-mean Gaussian variable with variance cr^ satisfies the bounds (16) and (17). In addition to the 
Gaussian case, the class of sub-Gaussian variates includes any bounded random variable (e.g., Bernoulli, multino- 
mial, uniform), any random variable with strictly log-concave density [6, 17], and any finite mixture of sub-Gaussian 
variables. 

The following lemma, proved in Appendix D, shows that the entries of the sample covariance based on i.i.d. 
samples of sub-Gaussian random vector satisfy an exponential-type tail bound with exponent a = 2. The argument 
is along the lines of a result due to Bickel and Levina [ I ], but with more explicit control of the constants in the error 
exponent: 

Lemma 1. Consider a zero-mean random vector [Xi , . . . , Xp) with covariance S* such that each Xi / -^/S^ is sub- 
Gaussian with parameter a. Given n i.i.d. samples, the associated sample covariance E" satisfies the tail bound 

F[|E^.-E*|>.] < 4exp{- ^^^^^^^^,^,^^^^^^^^^J , 
for all Se (0,max,;(S*)8(l + 4ct2)). 
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Thus, the sample covariance entries the tail condition T{f,Vt) with = [maxi(E*j) 8(1 + 4(t^)] ^, and an 
exponential-type tail function with a = 2 — namely 

f{n,6) = \exp{c^n6^), with = [l28(l + 4cr2)^ max(S*j2] (18) 

A little calculation shows that the associated inverse functions take the form 



Sf ir;n) = .hM^ and Mr;S) = (19) 
V c* n c,o^ 



2.3.2 Tail bounds with moment bounds 



In the following lemma, proved in Appendix E, we show that given i.i.d. observations from random variables with 
bounded moments, the sample covariance entries satisfy a polynomial-type tail bound. See the papers [26, for 
related results on tail bounds for variables with bounded moments. 



Lemma 2. Suppose there exists a positive integer m and scalar Km G K such that for i — 1, 



E 



(20) 



For i.i.d. samples {X^^^^"^^-^, the sample covariance matrix satisfies the bound 



>5\ < 



{m2"+i22'"(max, S* )2™ (K„, + 1)} 



(21) 



Thus, in this case, the sample covariance satisfies the tail condition T(/, v^,) with — 0, so that the bound holds 
for all S e (0, oo), and with the polynomial-type tail function 



f{n, S) = c,n'"(52™ where c, ^ l/{m2™+i22™(max, E* )2™ (A',„ + 1)}. 
Finally, a little calculation shows that in this case, the inverse tail functions take the form 

(r/c*)!/" 



(r/c*)i/2m 
Sf(n,r) = — , and nf{S,r) 



(52 



(22) 



(23) 



3 Main results and some consequences 

In this section, we state our main results, and discuss some of their consequences. We begin in Section 3.1 by 
stating some conditions on the true concentration matrix 0* required in our analysis, including a particular type of 
incoherence or irrepresentability condition. In Section 3.2, we state our first main result — namely. Theorem 1 on 
consistency of the estimator O, and the rate of decay of its error in elementwise norm. Section 3.3 is devoted 
to Theorem 2 on the model selection consistency of the estimator. Section 3.4 is devoted the relation between the 
log-determinant estimator and the ordinary Lasso (neighborhood-based approach) as methods for graphical model 
selection; in addition, we illustrate our irrepresentability assumption for some simple graphs. Finally, in Section 3.5, 
we state and prove some corollaries of Theorem 1, regarding rates in Frobenius and operator norms. 



3.1 Conditions on covariance and Hessian 

Our results involve some quantities involving the Hessian of the log-determinant barrier (6), evaluated at the true 
concentration matrix 8*. Using standard results on matrix derivatives [^], it can be shown that this Hessian takes the 
form 



r* 



v|.9(e) 



e=e* 



(24) 
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where ^ denotes the Kronecker matrix product. By definition, F* is a x matrix indexed by vertex pairs, so 
that entry T^^ ^.^ corresponds to the second partial derivative gg, "gQ^, . evaluated at = O*. When X has 
multivariate Gaussian distribution, then F* is the Fisher information of the model, and by standard results on cumulant 
functions in exponential families [5], we have the more specific expression F*^ ^.^ = cov{XjXk, X^Xm}. For 
this reason, F* can be viewed as an edge-based counterpart to the usual covariance matrix S*. 
We define the set of non-zero off-diagonal entries in the model concentration matrix Q*: 

E{e*) := {{t,3)eVxV \ i^j,e>:^^0}, (25) 

and let 5(6*) = {E{Q*) U {(1, 1), . . . , (p,p)} be the augmented set including the diagonal. We let 5'=(6*) denote 
the complement of S{<d*) in the set {1, ... ,p} x {1, . . . ,p}, corresponding to all pairs {£, m) for which Q},^ = 0. 
When it is clear from context, we shorten our notation for these sets to S and S"^, respectively. Finally, for any two 
subsets T and T' ofV x V, we use T^rp, to denote the |T| x \T'\ matrix with rows and columns of F* indexed by T 
and T' respectively. 

Our main results involve the ^oo/^oo norm applied to the covariance matrix S*, and to the inverse of a sub-block 
of the Hessian F*. In particular, we define 

p 

ifs* := iSlioo = (^.max (26) 

'■■■'^j=i 

corresponding to the £oo-operator norm of the true covariance matrix S*, and 

Kr^ msr'loc = |||([e*-i®e*-i]55)-^|||oo. (27) 

Our analysis keeps explicit track of these quantities, so that they can scale in a non-trivial manner with the problem 
dimension p. 

We assume the Hessian satisfies the following type of mutual incoherence or irrepresentable condition: 
Assumption 1. There exists some a e (0, 1] such that 

ll|r5cs(rS5)"'llloo < (l-a). (28) 

The underlying intuition is that this assumption imposes control on the influence that the non-edge terms, indexed 
by 5^, can have on the edge-based terms, indexed by 5*. It is worth noting that a similar condition for the Lasso, with 
the covariance matrix E* taking the place of the matrix F* above, is necessary and sufficient for support recovery using 
the ordinary Lasso [19, 22, 23, 26]. See Section 3.4 for illustration of the form taken by Assumption 1 for specific 
graphical models. 

A remark on notation: although our analysis allows the quantities K^.- , Ky- as well as the model size p and 
maximum node-degree d to grow with the sample size n, we suppress this dependence on n in their notation. 

3.2 Rates in elementwise £oo-norm 

We begin with a result that provides sufficient conditions on the sample size n for bounds in the elementwise l^o- 
norm. This result is stated in terms of the tail function /, and its inverses fif and 5f (equations (13) and (14)), and so 
covers a general range of possible tail behaviors. So as to make it more concrete, we follow the general statement with 
corollaries for the special cases of exponential-type and polynomial-type tail functions, corresponding to sub-Gaussian 
and moment-bounded variables respectively. 

In the theorem statement, the choice of regularization constant A„ is specified in terms of a user-defined parameter 
r > 2. Larger choices of r yield faster rates of convergence in the probability with which the claims hold, but also 
lead to more stringent requirements on the sample size. 
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Theorem 1. Consider a distribution satisfying the incoherence assumption (28) with parameter a G (0, 1], and the 
tail condition (12) with parameters T(f, f,). Let & be the unique optimum of the log-determinant program (11) with 
regularization parameter A„ = {8/a) Sf(n,p'^) for some t > 2. Then, if the sample size is lower bounded as 



(29) 



then with probability greater than 1 — l/p"^ ^ — > 1, we have: 

(a) The estimate Q satisfies the elementwise £oc-bound: 

||e-e*||oo < {2(l + 8«-i)A'r.} (30) 

(b) It specifies an edge set E{&) that is a subset of the true edge set E{Q*), and includes all edges with 
|e*,| > {2(l + 8a-i)Ar.} Sfin,p-). 

If we assume that the various quantities Kr',K^*,a remain constant as a function of {n,p,d), we have the 
elementwise too bound ||6 — 0*||oo = 0{Sf{n,p'^)), so that the inverse tail function Sf{n,p'^) (see equation (14)) 
specifies rate of convergence in the element-wise £oo-norm. In the following section, we derive the consequences of 
this ^oo -bound for two specific tail functions, namely those of exponential-type with a = 2, and polynomial-type tails 
(see Section 2.3). Turning to the other factors involved in the theorem statement, the quantities K^- and Kr* measure 
the sizes of the entries in the covariance matrix E* and inverse Hessian (F*)^^ respectively. Finally, the factor (1 + ^) 
depends on the irrepresentability assumption 1, growing in particular as the incoherence parameter a approaches 0. 

3.2.1 Exponential-type tails 

We now discuss the consequences of Theorem 1 for distributions in which the sample covariance satisfies an exponential- 
type tail bound with exponent a = 2. In particular, recall from Lemma 1 that such a tail bound holds when the variables 
are sub-Gaussian. 



Corollary 1. Under the same conditions as Theorem 1, suppose moreover that the variables Xi / ■^/S*^ are sub- 
Gaussian with parameter a, and the samples are drawn independently. Then if the sample size n satisfies the bound 

n > Ci d2(l + -)Mrlogp + log4) (31) 
a 

where Ci := {48^2 (1 + 4cr2) maxi(I]*j) max{K^tKr*, K^^Ky*}} , then with probability greater than 1 — 
the estimate Q satisfies the bound. 



|e-e*||oo < {l6\/2(l + 4a2) max(S*)(l + 8a-i)A>} J^^^^^i^ 



Proof. From Lemma 1, when the rescaled variables Xi/ ^JTjf^ are sub-Gaussian with parameter cr, the sample covari- 
ance entries satisfies a tail bound T(/, v^, ) with with = [ maxj 

where c* = [l28(l + 4(t^)^ maxi(S*j)^] ^ . As a consequence, for this particular model, the inverse functions 5/ (rtjp"^) 
and nf{6,p'^) take the form 



Sfin^pn = J'^^^ - ■/l28(l + 4.^)^max(S-)^./ -^°g^ + ^°g^ (32a) 
y n V « V n 

log(4p^) ^ ^ 2n2 .^*^2 /''rl0gP + l0g4 



nfiS^p^) = ' ^ U8{l + Aaym^x{j:*,y ""^^ " j. (32b) 

Substituting these forms into the claim of Theorem 1 and doing some simple algebra yields the stated corollary. □ 
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When Kp-'jK^tja remain constant as a function of {n,p,d), the corollary can be summarized succinctly as 

a sample size of n = logp) samples ensures that an elementwise £oc bound \\Q — 6*||oo = ^i\J^^^) holds 
with high probability. In practice, one frequently considers graphs with maximum node degrees d that either remain 
bounded, or that grow sub-linearly with the graph size (i.e., d = o{p)). In such cases, the sample size allowed by the 
corollary can be substantially smaller than the graph size, so that for sub-Gaussian random variables, the method can 
succeed in the p ^ n regime. 

3.2.2 Polynomial-type tails 

We now state a corollary for the case of a polynomial-type tail function, such as those ensured by the case of random 
variables with appropriately bounded moments. 

Corollary 2. Under the assumptions of Theorem 1, suppose the reseated variables Xi/ y/T,*^ have 4m*'' moments 
upper bounded by Km, and the sampling is i.i.d. Then if the sample size n satisfies the bound 

n > C2rf'(l + -)'p"/™, (33) 

where C2 '■= {l2m[m{Km + 1)]'^ maxi{Y,*^) ma.x{K^tKY* , K^tK^,}^^, then with probability greater than 
1 — the estimate Q satisfies the bound, 

lie-e*||oo < {4m[m(A™ + l)]5k (l + i)A>} 

a 

Proof. Recall from Lemma 2 that when the rescaled variables Xi / -^E*- have bounded 4771*'' moments, then the 
sample covariance E satisfies the tail condition T(/, v*) with = 0, and with f(n, 5) = Ct:n™S'^™ with c* defined 
as c* = l/{m^™+^2'^'"(maxi E*;)^™ {Km + 1)}. As a consequence, for this particular model, the inverse functions 
take the form 

{2m[m{Km + 1)]^ maxE* } \ ^ , (34a) 

i V n 

{2m[m{Km + l)]^ maxE* (^). (34b) 

The claim then follows by substituting these expressions into Theorem 1 and performing some algebra. □ 

When the quantities {Kr* , K-^^. , a) remain constant as a function of {n,p,d). Corollary 2 can be summarized 
succinctly as n = 

^(^(fpT/m^ samples are sufficient to achieve a convergence rate in elementwise -norm of the 

order IjO — O* ||oc = ^(■\/~^)' with high probability. Consequently, both the required sample size and the rate of 
convergence of the estimator are polynomial in the number of variables p. It is worth contrasting these rates with the 
case of sub-Gaussian random variables, where the rates have only logarithmic dependence on the problem size p. 

3.3 Model selection consistency 

Part (b) of Theorem 1 asserts that the edge set E{Q) returned by the estimator is contained within the true edge set 
E{Q*) — meaning that it correctly excludes all non-edges — and that it includes all edges that are "large", relative to 
the 5f{n,p'^) decay of the error. The following result, essentially a minor refinement of Theorem 1 , provides sufficient 
conditions linking the sample size n and the minimum value 

Oram mill 1 9* I (35) 

for model selection consistency. More precisely, define the event 

X(e;e*) := {sign(e,,) = sign(e*,.) V(z,j)ei?(e*)} (36) 




5f{n,p'') 
nf{5,p"') 



\l/2m 



\ 1/m 



<52 
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that the estimator 9 has the same edge set as 8*, and moreover recovers the correct signs on these edges. With this 
notation, we have: 

Theorem 2. Under the same conditions as Theorem 1, suppose that the sample size satisfies the lower bound 

n > n/^l/max{2A>.(l + 8a^^)6l-^[„, i;,, 6(l + 8a-i) dmax{A's-ifr*,i^|-i^r4 (37) 

Then the estimator is model selection consistent with high probability as p oo, 

f[M{%]Q*)\ > ^ i_ (38) 

2K * (1+ — ) 

In comparison to Theorem 1, the sample size requirement (37) differs only in the additional term — ^-^-^ — — 
involving the minimum value. This term can be viewed as constraining how quickly the minimum can decay as a 
function of {n,p), as we illustrate with some concrete tail functions. 



3.3.1 Exponential-type tails 

Recall the setting of Section 2.3.1, where the random variables {X^'-'^^/^I]*-} are sub-Gaussian with parameter a. 
Let us suppose that the parameters {Kr* , K^,* , a) are viewed as constants (not scaling with {p, d). Then, using the 
expression (32) for the inverse function in this setting, a corollary of Theorem 2 is that a sample size 

n - n{{d' + 9-^l)T\ogp) (39) 

is sufficient for model selection consistency with probability greater than 1 — Alternatively, we can state that 

n = Q{Td'^ logp) samples are sufficient, as along as the minimum value scales as 6„^in = ^( V ^^^). 



3.3.2 Polynomial-type tails 

Recall the setting of Section 2.3.2, where the rescaled random variables X^/y^E*^ have bounded 4m*'' moments. 
Using the expression (34) for the inverse function fif in this setting, a corollary of Theorem 2 is that a sample size 

n = niid' + O-Dp^/"-) (40) 

is sufficient for model selection consistency with probability greater than 1 — 1 /p^^^. Alternatively, we can state than 

n = ^l{d^p'^/™') samples are sufficient, as long as the minimum value scales as 6^^^^^ ~ Vl{p'^ Z'^'^"^'' / ^/n). 



3.4 Comparison to neighbor-based graphical model selection 

Suppose that X follows a multivariate Gaussian distribution, so that the structure of the concentration matrix 8* spec- 
ifies the structure of a Gaussian graphical model. In this case, it is interesting to compare our sufficient conditions for 
graphical model consistency of the log-determinant approach, as specified in Theorem 2, to those of the neighborhood- 
based method, first proposed by Meinshausen and Biihlmann [19]. The latter method estimates the full graph structure 
by performing an -regularized linear regression (Lasso) — of the form Xi = X^j^^i ^ij-^j + ^ — of ediCh node on 
its neighbors and using the support of the estimated regression vector 9 to predict the neighborhood set. These neigh- 
borhoods are then combined, by either an OR rule or an AND rule, to estimate the full graph. Various aspects of the 
high-dimensional model selection consistency of the Lasso are now understood [19, 23, 26]; for instance, it is known 
that mutual incoherence or irrepresentability conditions are necessary and sufficient for its success [22, 26]. In terms 
of scaling, Wainwright [23] shows that the Lasso succeeds with high probability if and only if the sample size scales 
as n x c{{d + 6',7ifn} logp)' where c is a constant determined by the covariance matrix S*. By a union bound over 
the p nodes in the graph, it then follows that the neighbor-based graph selection method in turn succeeds with high 
probability if n = rj({d + 9^^^} \ogp). 
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For comparison, consider the application of Theorem 2 to the case where the variables are sub-Gaussian (which 
includes the Gaussian case). For this setting, we have seen that the scaling required by Theorem 2 is n = + 
^minl logp)' so that the dependence of the log-determinant approach in 9„^in is identical, but it depends quadratically 
on the maximum degree d. We suspect that that the quadratic dependence might be an artifact of our analysis, but 
have not yet been able to reduce it to d. Otherwise, the primary difference between the two methods is in the nature 
of the irrepresentability assumptions that are imposed: our method requires Assumption 1 on the Hessian F*, whereas 
the neighborhood-based method imposes this same type of condition on a set of p covariance matrices, each of size 
(p— 1) X (p— 1), one for each node of the graph. Below we show two cases where the Lasso irrepresentability condition 
holds, while the log-determinant requirement fails. However, in general, we do not know whether the log-determinant 
irrepresentability strictly dominates its analog for the Lasso. 

3.4.1 Illustration of irrepresentability: Diamond graph 

Consider the following Gaussian graphical model example from Meinshausen [18]. Figure 2(a) shows a diamond- 
shaped graph G = {V, E), with vertex set V = {1,2,3,4} and edge-set as the fully connected graph over V with 
the edge (1,4) removed. The covariance matrix E* is parameterized by the correlation parameter p e [0, 1/V2]: 



1 




3 4 
(a) (b) 



Figure 2: (a) Graph of the example discussed by Meinshausen [18]. (b) A simple 4-node star graph. 

the diagonal entries are set to E* — 1, for all i G V; the entries corresponding to edges are set to E*^ = p for 
G E\{{2, 3)}, E23 = 0; and finally the entry corresponding to the non-edge is set as E*4 = 2p^. Meinshausen 
[ 1 8] showed that the £1 -penalized log-determinant estimator 6 fails to recover the graph structure, for any sample size, 
if p > — 1 + (3/2)^/^ « 0.23. It is instructive to compare this necessary condition to the sufficient condition provided 
in our analysis, namely the incoherence Assumption 1 as applied to the Hessian T*. For this particular example, a 
little calculation shows that Assumption 1 is equivalent to the constraint 

4|p|(|p| + l) < 1, 

an inequality which holds for all p e (—0.2017, 0.2017). Note that the upper value 0.2017 is just below the necessary 
threshold discussed by Meinshausen [18]. On the other hand, the irrepresentability condition for the Lasso requires 
only that 2\p\ < 1, i.e., p G (—0.5, 0.5). Thus, in the regime \p\ € [0.2017, 0.5), the Lasso irrepresentability condition 
holds while the log-determinant counterpart fails. 

3.4.2 Illustration of irrepresentability: Star graphs 

A second interesting example is the star-shaped graphical model, illustrated in Figure 2(b), which consists of a single 
hub node connected to the rest of the spoke nodes. We consider a four node graph, with vertex set V — {1, 2, 3, 4} 
and edge-set E = {{l,s) \ s G {2,3,4}}. The covariance matrix E* is parameterized the correlation parameter 
p E [—1,1]: the diagonal entries are set to E*j = 1, for all i E V; the entries corresponding to edges are set to 
E* = p for {i, i) e E; while the non-edge entries are set as E* = p^ for (i, j) ^ E. Consequently, for this particular 
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example. Assumption 1 reduces to the constraint + 2) < 1, which holds for all p E (—0.414,0.414). The 

irrepresentability condition for the Lasso on the other hand allows the full range p E ( — 1,1). Thus there is again 
a regime, \p\ G [0.414, 1), where the Lasso irrepresentability condition holds while the log-determinant counterpart 
fails. 

3.5 Rates in Frobenius and spectral norm 

We now derive some corollaries of Theorem 1 concerning estimation of 6* in Frobenius norm, as well as the spectral 
norm. Recall that s = [^^(O* ) | denotes the total number of off-diagonal non-zeros in Q*. 

Corollary 3. Under the same assumptions as Theorem 1, with probability at least 1 — the estimator 8 satisfies 

|||e-e*|||F < {2A'r.(l + ^)} VF+^(5/(n,p^), and (41a) 
|||e-e*|||2 < {2is:r-(l + ^)} min{^/F+^, d}(S/(n,p^). (41b) 

Proof. With the shorthand notation ly := 2Kr*{l + 8/a) Sf{n,p'^), Theorem 1 guarantees that, with probability at 
least 1 — l/p'^~^, II 6 — 9*11 oo < i^- Since the edge set of 6 is a subset of that of 8*, and 6* has at most p + s 
non-zeros (including the diagonal), we conclude that 

ll|e-e*|||^^ = [^(e„;-e:,)2+ 

< V y/s +p, 

from which the bound (41a) follows. On the other hand, for a symmetric matrix, we have 

lie -6*1 2 < lie -e* loo < diy, (42) 

using the definition of the i^oo-operator norm, and the fact that e and Q* have at most d non-zeros per row. Since the 
Frobenius norm upper bounds the spectral norm, the bound (41b) follows. 

□ 



3.5.1 Exponential-type tails 

For the exponential tail function case where the rescaled random variables Xi j -^E*- are sub-Gaussian with parameter 
a, we can use the expression (32) for the inverse function 5/ to derive rates in Frobenius and spectral norms. When the 
quantities Ky* , K^- , a remain constant, these bounds can be summarized succinctly as a sample size n — il{d^ logp) 
is sufficient to guarantee the bounds 



|e-ei|b = 0[J-^±^^], and (43a) 



|e-eib = 0( >"H^+P- d^}iogp ^43^^ 



with probability at least 1 — 1/p 



3.5.2 Polynomial-type tails 



r-2 



Similarly, let us again consider the polynomial tail case, in which the rescaled variates Xi/ ^S,*- have bounded 4m*'' 
moments and the samples are drawn i.i.d. Using the expression (34) for the inverse function we can derive rates in 
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the Frobenius and spectral norms. When the quantities Kr* , K^-' , a are viewed as constant, we are guaranteed that a 
sample size n ~ sufficient to guarantee the bounds 



||e-e*^ = ),and (44a) 



with probability at least 1 — 1 /p^^^. 

3.6 Rates for the covariance matrix estimate 

Finally, we describe some bounds on the estimation of the covariance matrix S*. By Lemma 3, the estimated concen- 
tration matrix Q is positive definite, and hence can be inverted to obtain an estimate of the covariance matrix, which 

we denote as E := (6)~^. 

Corollary 4. Under the same assumptions as Theorem 1, with probability at least 1 — l/p"^^^, the following bounds 
hold. 

(a) The element-wise norm of the deviation (E — E*) satisfies the bound 

||E-E*||oo < C3,[Sf{n,p^)]+Cid[Sf{n,p^)f (45) 
where C3 = 2KI,Kt* (l + f ) and C4 = QK^.K^. (l + f ) ^ 

(b) The £2 operator-norm of the deviation (E — E*) satisfies the bound 

|||S-E*|||2 < C,d[Sfin,p^)]+Cid^[Sfin,p^)f. (46) 

The proof involves certain lemmata and derivations that are parts of the proofs of Theorems 1 and 2, so that we 
defer it to Section 4.5. 

4 Proofs of main result 

In this section, we work through the proofs of Theorems 1 and 2. We break down the proofs into a sequence of lemmas, 
with some of the more technical aspects deferred to appendices. 

Our proofs are based on a technique that we call a primal-dual witness method, used previously in analysis of the 
Lasso [ : ]. It involves following a specific sequence of steps to construct a pair (0, Z) of symmetric matrices that 
together satisfy the optimality conditions associated with the convex program (11) with high probability. Thus, when 
the constructive procedure succeeds, Q is equal to the unique solution 9 of the convex program (11), and Z is an 
optimal solution to its dual. In this way, the estimator 6 inherits from 6 various optimality properties in terms of its 
distance to the truth 0*, and its recovery of the signed sparsity pattern. To be clear, our procedure for constructing 
is not a practical algorithm for solving the log-determinant problem (11), but rather is used as a proof technique for 
certifying the behavior of the M-estimator (11). 

4.1 Primal-dual witness approach 

As outlined above, at the core of the primal-dual witness method are the standard convex optimality conditions that 
characterize the optimum of the convex program (11). For future reference, we note that the sub-differential of the 
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norm || • || l off evaluated at some Q consists the set of all symmetric matrices Z S M.p^p such that 

{0 ifi^j 
sign(e.y ) ifi^j and 6,, ^ (47) 
e[-l,+l] ifi T^jandOy =0. 

The following result is proved in Appendix A; 

Lemma 3. For any A„ > and sample covariance E with strictly positive diagonal, the ti-regularized log-determinant 
problem (11) has a unique solution Q >- Q characterized by 

+ XnZ ^ 0, (48) 

where Z is an element of the subdijferential d\\Q\\i,ofi. 

Based on this lemma, we construct the primal-dual witness solution (8, Z) as follows: 

(a) We determine the matrix 8 by solving the restricted log-determinant problem 

8 := arg min {((8, S)) - logdGt(8) + A„||8||i,off}. (49) 

Note that by construction, we have 8 ;^ 0, and moreover 85c = 0. 

(b) We choose Zs as a member of the sub-differential of the regularizer || • || i,off , evaluated at 8. 

(c) We set Z5C as 

Zs^ = -^{-Ssc + [e-i]sc}, (50) 

which ensures that constructed matrices (8, Z) satisfy the optimality condition (48). 

(d) We verify the strict dual feasibility condition 

\Zij\ < 1 forall e 

To clarify the nature of the construction, steps (a) through (c) suffice to obtain a pair (8, Z) that satisfy the optimality 
conditions (48), but do nof guarantee that Z is an element of sub-differential 9||8||i.off ■ By construction, specifically 
step (b) of the construction ensures that the entries Z in S* satisfy the sub-differential conditions, since Zs is a member 
of the sub-differential of 9|18s|ji off- The purpose of step (d), then, is to verify that the remaining elements of Z 
satisfy the necessary conditions to belong to the sub-differential. 

If the primal-dual witness construction succeeds, then it acts as a witness to the fact that the solution 8 to the 
restricted problem (49) is equivalent to the solution 8 to the original (unrestricted) problem (11). We exploit this fact 
in our proofs of Theorems 1 and 2 that build on this: we first show that the primal-dual witness technique succeeds with 
high-probability, from which we can conclude that the support of the optimal solution 8 is contained within the support 
of the true 8*. In addition, we exploit the characterization of 8 provided by the primal-dual witness construction to 
establish the elementwise loo bounds claimed in Theorem 1. Theorem 2 requires checking, in addition, that certain 
sign consistency conditions hold, for which we require lower bounds on the value of the minimum value 6'inin- 

In the analysis to follow, some additional notation is useful. We let W denote the "effective noise" in the sample 
covariance matrix E, namely 

W := %-{Q*)-^. (51) 

Second, we use A = 8—8* to measure the discrepancy between the primal witness matrix 8 and the truth 8*. Finally, 
recall the log-determinant barrier g from equation (6). We let i?(A) denote the difference of the gradient V_g(8) = 
8^^ from its first-order Taylor expansion around 8*. Using known results on the first and second derivatives of the 
log-determinant function (see p. 641 in Boyd and Vandenberghe [3]), this remainder takes the form 

i?(A) = 8^1 -8*"V8*"^A8*"\ (52) 
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4.2 Auxiliary results 



We begin with some auxiliary lemmata, required in the proofs of our main theorems. In Section 4.2.1, we provide 
sufficient conditions on the quantities W and R for the strict dual feasibility condition to hold. In Section 4.2.2, we 
control the remainder term i?(A) in terms of A, while in Section 4.2.3, we control A itself, providing elementwise 
£oo bounds on A. In Section 4.2.4, we show that under appropriate conditions on the minimum value d„^in, the bounds 
in the earlier lemmas guarantee that the sign consistency condition holds. All of the analysis in these sections is 
deterministic in nature. In Section 4.2.5, we turn to the probabilistic component of the analysis, providing control 
of the noise W in the sample covariance matrix. Finally, the proofs of Theorems 1 and 2 follows by using this 
probabilistic control of W and the stated conditions on the sample size to show that the deterministic conditions hold 
with high probability. 



4.2.1 Sufficient conditions for strict dual feasibility 

We begin by stating and proving a lemma that provides sufficient (deterministic) conditions for strict dual feasibility 
to hold, so that || j| < 1. 

Lemma 4 (Strict dual feasibility). Suppose that 

maxjllW^llo,, ||i?(A)!U} < (53) 
Then the matrix Zs<^ constructed in step (c) satisfies ||2^5<:||oo < 1, and therefore = 0. 

Proof. Using the definitions (51) and (52), we can re-write the stationary condition (48) in an alternative but equivalent 
form 

e*-^Ae*^'^ + W - R{A) + XnZ =0. (54) 

This is a linear-matrix equality, which can be re-written as an ordinary linear equation by "vectorizing" the matrices. 
We use the notation vcc{A), or equivalently A for the p^-vector version of the matrix A £ MP^p, obtained by stacking 
up the rows into a single column vector. In vectorized form, we have 

vec(e*"^Ae*"^) = (6*^^ ® e*"^)A = r*A. 

In terms of the disjoint decomposition S and 5^, equation (54) can be re-written as two blocks of linear equations as 
follows: 

T*ssAs + Ws-Rs + XnZs = (55a) 
r*s.s^s + Ws^-Rs^+KZs^- = 0. (55b) 

Here we have used the fact that A^e = by construction. 

Since Tgg is invertible, we can solve for A5 from equation (55a) as follows: 

As = {T*ssy'[-Ws + Rs-KZs]. 
Substituting this expression into equation (55b), we can solve for Zgc as follows: 

= -— Fgc^As + — — — T^go 

An An An 

= -^n^s{nsy\ws-Rs)+n.s{nsy'zs-^{ws^^Rs^)- (56) 

An An 
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Taking the £00 norm of both sides yields 

An 



l|rS=5(rss) 'lIUII^slloo + T^(||T^5||oo + ||i?5||oo). 

A 11 



Recalhng Assumption 1 — namely, that 111^5=5 (FJ.^) ^|||oo < (1 — ot) — we have 

- 2 — a - 

ll^s^lU < ^— (||1^5||oo + ||i?s||oo) + (l-a), 

where we have used the fact that j|Z5||oo < 1, since Z belongs to the sub-differential of the norm || • ||i,off by 
construction. Finally, applying assumption (53) from the lemma statement, we have 

117 II ^ (2-a)/aA„ 



< | + < 1, 



as claimed. 



□ 



4.2.2 Control of remainder term 

Our next step is to relate the behavior of the remainder term (52) to the deviation A = — O*. 

Lemma 5 (Control of remainder). Suppose that the elementwise i^o bound \ \A\\oo < 3K^,d ^'^^Ids. Then: 

i?(A) = e*"^Ae*~^Aje*"\ (57) 

where J := Er=o(-l)'' l©*"^^) 

norm ||| J^|||oo ^ 3/2. Moreover, in terms of the elementwise loo-norm, we 

have 

\\Rm\oo < ^dWAWlKl,,. (58) 

We provide the proof of this lemma in Appendix B. It is straightforward, based on standard matrix expansion 
techniques. 

4.2.3 Sufficient conditions for £^0 bounds 

Our next lemma provides control on the deviation A — Q — 0*, measured in elementwise £00 norm. 
Lemma 6 (Control of A). Suppose that 

Then we have the elementwise £00 bound 

||A||oo - ||e-e*||oo < r. (60) 

We prove the lemma in Appendix C; at a high level, the main steps involved are the following. We begin by noting 
that Ggc = = 0, so that ||A||oc = HA^Hoo- Next, we characterize 85 in terms of the zero-gradient condition 
associated with the restricted problem (49). We then define a continuous map F : As ^ F{As) such that its fixed 
points are equivalent to zeros of this gradient expression in terms of A5 — Qs ~ ©s- We then show that the function 
F maps the ^oo-ball 

B(r) ~ {05 I lieslloo <r}, withr:=2Xr.(||W^||oo + A„), (61) 

onto itself. Finally, with these results in place, we can apply Brouwer's fixed point theorem (e.g., p. 161; Ortega and 
Rheinboldt [20]) to conclude that F does indeed have a fixed point inside B(r). 
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4.2.4 Sufficient conditions for sign consistency 

We now show how a lower bound on the minimum value 6',nin, when combined with Lemma 6, allows us to guarantee 
sign consistency of the primal witness matrix Qs- 

Lemma 7 (Sign Consistency). Suppose the minimum absolute value O^^^ of non-zero entries in the true concentration 
matrix O* is lower bounded as 

0min > 4Xr.(||VF||oo + A„), (62) 

then sign(Qs) = sign(Qg) holds. 

This claim follows from the bound (62) combined with the bound (60) ,which together imply that for all G S, 
the estimate Qij cannot differ enough from Q*j to change sign. 

4.2.5 Control of noise term 

The final ingredient required for the proofs of Theorems 1 and 2 is control on the sampling noise = S — S*. This 
control is specified in terms of the decay function / from equation (12). 



Lemma 8 (Control of Sampling Noise). For any r > 2 and sample size n such that Sf{n,p'^) < 1 /v,, we have 

\W\\.^ > Sfin^p-") 



< ^ ^ 0. (63) 



Proof. Using the definition (12) of the decay function /, and applying the union bound over all p^ entries of the noise 
matrix, we obtain that for all (5 < l/v^. 



max I 



■-m,\>5\ < pVf{n,S). 
Setting S — Sf{n,p'^) yields that 



m^x\W,,\>Sfin,p-)] < py[f{njf{n,p-))] ^ l/p-~^ 



as claimed. Here the last equality follows since /(n, 6f{n,p'^)) = p'^ , using the definition (14) of the inverse function 
5f. □ 

4.3 Proof of Theorem 1 

We now have the necessary ingredients to prove Theorem 1. We first show that with high probability the witness 
matrix Q is equal to the solution 8 to the original log-determinant problem (11), in particular by showing that the 
primal-dual witness construction (described in in Section 4.1) succeeds with high probability. Let A denote the event 
that ||VK||oo < Sf{n,p'^). Using the monotonicity of the inverse tail function (15), the lower lower bound (29) on the 
sample size n implies that 6f{n,p'^) < Consequently, Lemma 8 implies that P(^) > 1 — Accordingly, 

we condition on the event A in the analysis to follow. 

We proceed by verifying that assumption (53) of Lemma 4 holds. RecalUng the choice of regularization penalty 
A„ = (8/a) 6f{n,p'^), we have IjM^lloo < (a/8)A„. In order to establish condition (53) it remains to establish the 
bound |li?(A)||oo < ^F^- '^'^ '^^^ steps, by using Lemmas 6 and 5 consecutively. First, we show that the 
precondition (59) required for Lemma 6 to hold is satisfied under the specified conditions on n and A„ . From Lemma 8 
and our choice of regularization constant A„ = (8/ a) 6f{n,p'^), 

2i^r-(||VP^||oo + A„) < 2Kr^(l + -)Sfin,p^), 
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provided 5f{n,p'^) < From the lower bound (29) and the monotonicity (15) of the tail inverse functions, we 

have 

2Kr^{l + l.)s,in,pn < -.{^,^^3;^}, (64) 
showing that the assumptions of Lemma 6 are satisfied. Applying this lemma, we conclude that 

||A|U < 2Xr.(!|W^||oo + A„) < 2Kr*(l + ^)Sfin,p^). (65) 

Turning next to Lemma 5, we see that its assumption || Aj|oo < 3k ,d holds, by applying equations (64) and (65). 
Consequently, we have 

||i?(A)|U < ld\\A\\lK% 

< 6K%K^,d{l + ^y[Sf{n,pnr 



as required, where the final inequality follows from our condition (29) on the sample size, and the monotonicity 
property (15). ^ 

Overall, we have shown that the assumption^(53) of Lemma 4 holds, allowing us to conclude that 9 = 6. The 
estimator Q then satisfies the £oo-bound (65) of O, as claimed in Theorem 1(a), and moreover, we have Qs" = ©s^ = 
0, as claimed in Theorem 1(b). Since the above was conditioned on the event A, these statements hold with probabihty 
P(^) > 1 



4.4 Proof of Theorem 2 

We now turn to the proof of Theorem 2. A little calculation shows that the assumed lower bound (37) on the sample 
size n and the monotonicity property (15) together guarantee that 

/ 8 \ - 

Proceeding as in the proof of Theorem 1, with probability at least 1 — we have the equality = 0, and 

also that ||0 — 0*||oo < 6*111111/2. Consequently, Lemma 7 can be applied, guaranteeing that sign(0*j) = sign(Oy ) 
for all G E. Overall, we conclude that with probability at least 1 — the sign consistency condition 

sign(0* ) = sign(0.y) holds for all G E, as claimed. 



4.5 Proof of Corollary 4 

With the shorthand A = O — O*, we have 

S - E* = (0* + A)-^ - {e*)~\ 

From the definition (52) of the residual this difference can be written as 

E-S* = -0*"^A0*"^ +i?(A). (66) 
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Proceeding as in the proof of Theorem 1 we condition on the event A = {HW^Hcxj < ^/("^iP^)}' which 
holds with probability P(^) > 1 — — As in the proof of that theorem, we are guaranteed that the assumptions of 
Lemma 5 are satisfied, allowing us to conclude 

i?(A) = e*"^Ae*"^Aje*"\ (67) 

where J := T,T=o(-^)'' i^*'^^)'' has norm ||| J^IU < 3/2. 

We begin by proving the bound (45). From equation (66), we have — E*|loo < |lI/(A)||oo + |li?(A)||oo- From 
Lemma 5, we have the elementwise £oo-norm bound 

||i?(A)||oo < ld\\A\\lKl,. 

The quantity L{A) in turn can be bounded as follows, 

||L(A)||oo = max|efe*"^Ae*"^ej| 



< max||efe*"^||imax||Ae*"^ 



3 



< max||efe*"'l|il|A||oo||max|ie*"'e,||i 



where we used the inequality that ||Au|joo < ||A||oo||u||i- Simplifying further, we obtain 

||L(A)|u < |||e*-^||U||A||^|||e*-i|||i 

< |||e*-^f^||A|u 

< ^i.||A||oo, 

where we have used the fact that |||6*~^|||i = III [Q*^^]^ll|oo = ||| 6* ~^|||oo> which follows from the symmetry of 9*"^. 
Combining the pieces, we obtain 

||S-I]*i|oo < i|L(A)||oo + ||i?(A)||oo (68) 

The claim then follows from the elementwise £oo-norm bound (30) from Theorem I. 

Next, we establish the bound (46) in spectral norm. Taking the £oo operator norm of both sides of equation (66) 

yields the inequahty |||I] — 2]*|||oo < |||^(A)|||oo + |||-R(A)|||oo- Using the expansion (67), and the sub-multiplicativity 
of the £oo operator norm, we obtain 



ll|i?(A)|||oo < |||e*-i|||oo|||A||uil|e*-i|||c 
11^ 

3 



oo <^ oo 



< llie*-^|||LllUllloolllA|||2 



2 



where the last inequality uses the bound ||| J|||oo < 3/2. (Proceeding as in the proof of Lemma 5, this bound holds 
conditioned on A, and for the sample size specified in the theorem statement.) In turn, the term L( A) can be bounded 

as 

|||i(A)||u < |||e*-iAe*-^||U 

< |||e*-^|||Lll|A|||oo 

< if|.|||A||U, 
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(a) (b) (c) 

Figure 3. Illustrations of different graph classes used in simulations, (a) Chain (rf = 2). (b) Four-nearest neighbor grid 
(d = 4) and (c) Star-shaped graph (d G {1, . . . , p — 1}). 

where the second inequahty uses the sub-muhipUcativity of the £oo -operator norm. Combining the pieces yields 

|||£-I]*||U < |||L(A)||U + |||i?(A)||U < At.|||A||U + ^if|.||A||^ (69) 
Conditioned on the event A, we have the bound (42) on the -operator norm 

|||A||U < 2Kr,(l + ^)d5f{n,p'). 

Substituting this bound, as well as the elementwise £oo-norm bound (30) from Theorem 1, into the bound (69) yields 
the stated claim. 

5 Experiments 

In this section, we illustrate our results with various experimental simulations, reporting results in terms of the proba- 
bility of correct model selection (Theorem 2) or the ^oo-error (Theorem 1). For these illustrations, we study the case 
of Gaussian graphical models, and results for three different classes of graphs, namely chains, grids, and star-shaped 
graphs. We also consider various scalings of the quantities which affect the performance of the estimator: in addition 
the triple {n, p, d), we also report some results concerning the role of the parameters Ks* , Kr* and d^in that we have 
identified in the main theorems. For all results reported here, we solved the resulting ^i-penalized log-determinant 
program (11) using the glasso program of Friedman et al. [10], which builds on the block co-ordinate descent 
algorithm of d' Aspremont et al. [8]. 

Figure 3 illustrates the three types of graphs used in our simulations: chain graphs (panel (a)), four-nearest neighbor 
lattices or grids (panel (b)), and star-shaped graphs (panel (c)). For the chain and grid graphs, the maximal node degree 
d is fixed by definition, to d = 2 for chains, and d = 4 for the grids. Consequently, these graphs can capture the 
dependence of the required sample size n only as a function of the graph size p, and the parameters (if s. , Kr* , d„^in)- 
The star graph allows us to vary both d and p, since the degree of the central hub can be varied between 1 and p — 1. 
For each graph type, we varied the size of the graph p in different ranges, from p = 64 upwards to p = 375. 

For the chain and star graphs, we define a covariance matrix S* with entries S]*j = 1 for all z = 1, . . . ,p, and 
= p for all € E for specific values of p specified below. Note that these covariance matrices are sufficient to 
specify the full model. For the four-nearest neighbor grid graph, we set the entries of the concentration matrix O*^ = oj 
for (i, j) G E, with the value uj specified below. In all cases, we set the regularization parameter A„ proportional to 
^y\og{p)/n, as suggested by Theorems 1 and 2, which is reasonable since the main purpose of these simulations 
is to illustrate our theoretical results. However, for general data sets, the relevant theoretical parameters cannot be 
computed (since the true model is unknown), so that a data-driven approach such as cross-validation might be required 
for selecting the regularization parameter A„. 
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Given a Gaussian graphical model instance, and the number of samples n, we drew = 100 batches of n 
independent samples from the associated multivariate Gaussian distribution. We estimated the probability of correct 
model selection as the fraction of the TV = 100 trials in which the estimator recovers the signed-edge set exactly. 

Note that any multivariate Gaussian random vector is sub-Gaussian; in particular, the rescaled variates Xi/ -^E*- 
are sub-Gaussian with parameter cr = 1, so that the elementwise ^oo-bound from Corollary 1 applies. Suppose we 
collect relevant parameters such as ^min and the covariance and Hessian-related terms Ks* , Kr* and a into a single 
"model-complexity" term K defined as 



K := 

Then, as a corollary of Theorem 2, a sample size of order 



(1 + 8a" ^) (max S* ) max{ifs*ifr* , if|.i^r* , 



(70) 



n = n{K^ (frlogp) , (71) 

is sufficient for model selection consistency with probability greater than 1 — 1/p^"^. In the subsections to follow, we 
investigate how the empirical sample size n required for model selection consistency scales in terms of graph size p, 
maximum degree d, as well as the "model-complexity" term K defined above. 




Figure 4. Simulations for chain graphs with varying number of nodes p, edge covariances S*j = 0.10. Plots of probability 
of correct signed edge-set recovery plotted versus the ordinary sample size n in panel (a), and versus the rescaled sample 
size n/ logp in panel (b). Each point corresponds to the average over 100 trials. 



5.1 Dependence on graph size 

Panel (a) of Figure 4 plots the probability of correct signed edge-set recovery against the sample size n for a chain- 
structured graph of three different sizes. For these chain graphs, regardless of the number of nodes p, the maximum 
node degree is constant d = 2, while the edge covariances are set as = 0.2 for all {i, j) E E, so that the quantities 
{Kt,* , Kr* , a) remain constant. Each of the curve in panel (a) corresponds to a different graph size p. For each curve, 
the probability of success starts at zero (for small sample sizes n), but then transitions to one as the sample size is 
increased. As would be expected, it is more difficult to perform model selection for larger graph sizes, so that (for 
instance) the curve for p = 375 is shifted to the right relative to the curve for p = 64. Panel (b) of Figure 4 replots 
the same data, with the horizontal axis rescaled by (1/ logp). This scaling was chosen because for sub-Gaussian tails, 
our theory predicts that the sample size should scale logarithmically with p (see equation (71)). Consistent with this 
prediction, when plotted against the rescaled sample size n/ \ogp, the curves in panel (b) all stack up. Consequently, 
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Star graph Star graph 




40 



Figure 5. Simulations for a star graph with varying number of nodes p, fixed maximal degree d — 40, and edge covariances 
T,*j — 1/16 for all edges. Plots of probability of correct signed edge-set recovery versus the sample size n in panel (a), 
and versus the reseated sample size n/ log p in panel (b). Each point corresponds to the average over A'^ — 100 trials. 



the ratio {n/ logp) acts as an effective sample size in controlling the success of model selection, consistent with the 
predictions of Theorem 2 for sub-Gaussian variables. 

Figure 5 shows the same types of plots for a star-shaped graph with fixed maximum node degree d = 40, and 
Figure 6 shows the analogous plots for a grid graph with fixed degree d = 4. As in the chain case, these plots show 
the same type of stacking effect in terms of the scaled sample size n/ logp, when the degree d and other parameters 
((a, A'r* , K^* )) are held fixed. 

5.2 Dependence on the maximum node degree 

Panel (a) of Figure 7 plots the probability of correct signed edge-set recovery against the sample size n for star- 
shaped graphs; each curve corresponds to a different choice of maximum node degree d, allowing us to investigate the 
dependence of the sample size on this parameter. So as to control these comparisons, the models are chosen such that 
quantities other than the maximum node-degree d are fixed: in particular, we fix the number of nodes p = 200, and the 
edge covariance entries are set as 1]*^ ~ 2.5/ d for {i, j) G E so that the quantities {K^* , Kr* , a) remain constant. 
The minimum value 6'min in turn scales as 1 /d. Observe how the plots in panel (a) shift to the right as the maximum 
node degree d is increased, showing that star-shaped graphs with higher degrees are more difficult. In panel (b) of 
Figure 7, we plot the same data versus the rescaled sample size n/d. Recall that if all the curves were to stack up 
under this rescaling, then it means the required sample size n scales linearly with d. These plots are closer to aligning 
than the unrescaled plots, but the agreement is not perfect. In particular, observe that the curve d (right-most in panel 
(a)) remains a bit to the right in panel (b), which suggests that a somewhat more aggressive rescaling — ^perhaps n/d"' 
for some 7 G (1, 2) — is appropriate. 

Note that for 6',nin scaling as 1/d, the sufficient condition from Theorem 2, as summarized in equation (71), is 
n = ri((i^ logp), which appears to be overly conservative based on these data. Thus, it might be possible to tighten 
our theory under certain regimes. 

5.3 Dependence on covariance and Hessian terms 

Next, we study the dependence of the sample size required for model selection consistency on the model complexity 
term K defined in (70), which is a collection of the quantities A's*, Kr* and a defined by the covariance matrix and 
Hessian, as well as the minimum value 0min- Figure 8 plots the probability of correct signed edge-set recovery versus 
the sample size n for chain graphs. Here each curve corresponds to a different setting of the model complexity factor 
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Figure 6. Simulations for 2-dimensional lattice with 4-nearest-neighbor interaction, edge strength interactions 0*j = 0.1, 
and a varying number of nodes p. Plots of probability of correct signed edge-set recovery versus the sample size n in panel 
(a), and versus the rescaled sample size n/ log p in panel (b). Each point corresponds to the average over TV = 100 trials. 



Truncated Star with Varying d Truncated Star with Varying d 




n/d 



(a) (b) 

Figure 7. Simulations for star graphs with fixed number of nodes p = 200, varying maximal (hub) degree d, edge 
covariances E*j = 2.5/d. Plots of probability of correct signed edge-set recovery versus the sample size n in panel (a), 
and versus the rescaled sample size n/d in panel (b). 
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Chain graph with varying K 




-»-K=64.1 
— K=44.9 
-e-K=37.5 



K=30.2 
— K=26.7 
K=22.9 



500 



n 



Figure 8. Simulations for cliain grapli witli fixed number of nodes p = 120, and varying model complexity K. Plot of 
probability of correct signed edge-set recovery versus the sample size n. 

K, but with a fixed number of nodes p = 120, and maximum node-degree d = 2. We varied the actor K by varying 
the value p of the edge covariances = p, G E. Notice how the curves, each of which corresponds to a 

different model complexity factor, shift rightwards as K is increased so that models with larger values of K require 
greater number of samples n to achieve the same probability of correct model selection. These rightward-shifts are in 
qualitative agreement with the prediction of Theorem 1, but we suspect that our analysis is not sharp enough to make 
accurate quantitative predictions regarding this scaling. 

5.4 Convergence rates in elementwise ^oo-norm 

Finally, we report some simulation results on the convergence rate in elementwise ^oo-norm. According to Corollary 1, 



in the case of sub-Gaussian tails, if the elementwise ^oo-norm should decay at rate C( y -^^)' once the sample size 
n is sufficiently large. Figure 9 shows the behavior of the elementwise £oo-norm for star-shaped graphs of varying 
sizes p. The results reported here correspond to the maximum degree d = \0.lp]; we also performed analogous 
experiments for d = 0{\ogp) and d ~ 0{1), and observed qualitatively similar behavior The edge correlations 
were set as S*j — 2.5/(iforall {i,j) G i? so that the quantities (Ks* , A'r* , a) remain constant. With these settings, 
each curve in Figure 9 corresponds to a different problem size, and plots the elementwise -error versus the rescaled 
sample size n/ logp, so that we expect to see curves of the form f{t) = l/v^. The curves show that when the rescaled 
sample size (n/ logp) is larger than some threshold (roughly 40 in the plots shown), the elementwise £oo norm decays 



6 Discussion 

The focus of this paper is the analysis of the high-dimensional scaling of the -regularized log determinant prob- 
lem (11) as an estimator of the concentration matrix of a random vector Our main contributions were to derive 
sufficient conditions for its model selection consistency as well as convergence rates in both elementwise £oo-norm, as 
well as Frobenius and spectral norms. Our results allow for a range of tail behavior, ranging from the exponential-type 
decay provided by Gaussian random vectors (and sub-Gaussian more generally), to polynomial-type decay guaranteed 
by moment conditions. In the Gaussian case, our results have natural interpretations in terms of Gaussian Markov 
random fields. 

Our main results relate the i.i.d. sample size n to various parameters of the problem required to achieve consistency. 
In addition to the dependence on matrix size p, number of edges s and graph degree d, our analysis also illustrates the 




at the rate W -2sp^ which is consistent with Corollary 1. 
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Figure 9. Simulations for a star graphi withi varying number of nodes p, maximum node degree d = [O.lp] , edge covari- 
ances E*j — 2.5/d. Plot of the element-wise loo norm of the concentration matrix estimate error ||B — 0*||oo versus the 
rescaled sample size n/ log(p). 

role of other quantities, related to the structure of the covariance matrix S* and the Hessian of the objective function, 
that have an influence on consistency rates. Our main assumption is an irrepresentability or mutual incoherence 
condition, similar to that required for model selection consistency of the Lasso, but involving the Hessian of the log- 
determinant objective function (11), evaluated at the true model. When the distribution of X is multivariate Gaussian, 
this Hessian is the Fisher information matrix of the model, and thus can be viewed as an edge-based counterpart to 
the usual node-based covariance matrix We report some examples where irrepresentability condition for the Lasso 
hold and the log-determinant condition fails, but we do not know in general if one requirement dominates the other. In 
addition to these theoretical results, we provided a number of simulation studies showing how the sample size required 
for consistency scales with problem size, node degrees, and the other complexity parameters identified in our analysis. 

There are various interesting questions and possible extensions to this paper. First, in the current paper, we have 
only derived sufficient conditions for model selection consistency. As in past work on the Lasso [23], it would also 
be interesting to derive a converse result — namely, to prove that if the sample size n is smaller than some function of 
{p, d, s) and other complexity parameters, then regardless of the choice of regularization constant, the log-determinant 
method fails to recover the correct graph structure. Second, while this paper studies the problem of estimating a fixed 
graph or concentration matrix, a natural extension would allow the graph to vary over time, a problem setting which 
includes the case where the observations are dependent. For instance, Zhou et al. [27] study the estimation of the 
covariance matrix of a Gaussian distribution in a time-varying setting, and it would be interesting to extend results of 
this paper to this more general setting. 

Acknowledgements 

We thank Shuheng Zhou for helpful comments on an earlier draft of this work. Work was partially supported by NSF 
grant DMS -0605 165. Yu also acknowledges support from ARO W91 lNF-05-1-0104, NSFC-60628102, and a grant 
from MSRA. 

A Proof of Lemma 3 

In this appendix, we show that the regularized log-determinant program (11) has a unique solution whenever A„ > 0, 
and the diagonal of the sample covariance S" is strictly positive. By the strict convexity of the log-determinant 
barrier [3], if the minimum is attained, then it is unique, so that it remains to show that the minimum is achieved. If 
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A„ > 0, then by Lagrangian duality, the problem can be written in an equivalent constrained form: 

min {((6, S")) -logdct(e)| (72) 

for some C(A„) < +oo. Since the off-diagonal elements remain bounded within the ^i-ball, the only possible issue 
is the behavior of the objective function for sequences with possibly unbounded diagonal entries. Since any Q in 
the constraint set is positive-definite, its diagonal entries are positive, and hence bounded below by zero. Further, by 
Hadamard's inequality for positive definite matrices [ 1 2], we have log det 8 < X]r=i log Qu, so that 

f]e,,sr,- log dote > f]{e,,sr,-ioge,,;}. 

1=1 1=1 

As long as I]"^ > for each i = 1, . . . ,p, this function is coercive, meaning that it diverges to infinity for any sequence 
11(0*]^, . . . , 0pp)|| +00. Consequently, the minimum is attained. 

Returning to the penalized form (11), by standard optimality conditions for convex programs, a matrix S >- is 
optimal if and only belongs to the sub-differential of the objective, or equivalently if and only if there exists a matrix 
Z in the sub-differential of the off-diagonal norm || • || i,off such that 

as claimed. 



B Proof of Lemma 5 

We write the remainder in the form 

R{A) = {e* + A)-^ -e*^^ + e*^^Ae*'\ 

By sub-multiplicativity of the ||| • |||oo matrix norm, for any two pxp matrices A, B, we have |||^ -B|||oo < 
so that 

l||e*-'A||u < |||e*"^||UII|A||u 

< Ks' d\\A\\oo < 1/3, (73) 

where we have used the definition of A's* , the fact that A has at most d non-zeros per row/column, and our assumption 
||Aj|oo < l/(3^s*)- Consequently, we have the convergent matrix expansion 

(e* + A)"^ = (e*(/ + e*"^A))"^ 
^ {i + e*-'Ay\e*)'' 



k=0 



= e*"^ - e*"^Ae*"V^(-i)'=(e*"^A) (e*) 

fc=2 

= e*"^ - e*"^Ae*"^ + e*"^Ae*"^Aje*"\ 

whereJ = Er=o(-l)'(0*"'A)'. 

We now prove the bound (58) on the remainder as follows. Let denote the unit vector with 1 in position i and 
zeroes elsewhere. From equation (57), we have 

||i?(A)|U - max|efe*"^A e*"^Aje*"^ej| 

< max lief e*"^A||oo max IIG*"^ Aje*"^e JN, 
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which follows from the fact that for any vectors a,b E W, {a^bl < llallooll^H i- This in turn can be simplified as, 

||i?(A)||oo < max||efe*"'||i ||A||oo max||e*"'Aje*"'e,||i 
« j 

since for any vector u E R^, ||m^A||oo < i|| Aj|oo, where || Aj|oc is the elementwise i'oc-norm. Continuing on, we 
have 

||i?(A)||oo < |||e*-'|||oo ||A||oo|| |||e*"'Aje*-'|||i, 

where |||^|||i := max||j.||j^x is the £i -operator norm. Since |||j4|||i ~ |||A-^|||oo, we have 

||i?(A)|u < ||Ai|oo|||e*-'||U |||e*-V^Ae*-'||U (74) 

< ||A|Ui^s.|l|0*"'lllLlll-/^llloo|||A||U 
Recall that J = J^'kLoi^^)'' (6*^^ A)'^. By sub-multiplicativity of ||| • matrix norm, we have 

oo 

iiijT||| < V iiiAe*~^iii'= < - < - 

'"°"-fro i-llie-^||UII|A||u - 2' 

since |||0*~^|||oo|||A|||oo < 1/3 from equation (73). Substituting this in (74), we obtain 

l|i?(A)||oo < lj\A\\^ K^4e*''\iUm^ 

< ld\\A\\lKl„ 

where the final line follows since ||| A|||oo < A|joo, and since A has at most d non-zeroes per row/column. 



C Proof of Lemma 6 

By following the same argument as in Appendix A, we conclude that the restricted problem (49) has a unique optimum 
Q. Let Z be any member of the sub-differential of j| • j| i.ofr, evaluated at 6. By Lagrangian theory, the witness must 
be an optimum of the associated Lagrangian problem 

min {((9, S)) - logdct(e) + A„((e, Z))} . 

In fact, since this Lagrangian is strictly convex, Q is the only optimum of this problem. Since the log-determinant 
barrier diverges as Q approaches the boundary of the positive semi-definite cone, we must have ;^ 0. If we take 
partial derivatives of the Lagrangian with respect to the unconstrained elements 65, these partial derivatives must 
vanish at the optimum, meaning that we have the zero-gradient condition 

G(es) = -65^ + Ss + XnZs = 0. (75) 

To be clear, Q is the pxp matrix with entries in S equal to Qs and entries in S**^' equal to zero. Since this zero-gradient 
condition is necessary and sufficient for an optimum of the Lagrangian problem, it has a unique solution (namely, 85). 

Our goal is to bound the deviation of this solution from 85, or equivalently to bound the deviation A = 8 — 8*. 
Our strategy is to show the existence of a solution A to the zero-gradient condition (75) that is contained inside the 
ball B(r) defined in equation (61). By uniqueness of the optimal solution, we can thus conclude that — 8* belongs 
this ball. In terms of the vector A5 — 85 — 8*5, let us define a map F : RI'^I — > RI"^! via 

F(As) := -{T*ssy\G{Q*s + As))+Es, (76) 
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where G denotes the vectorized form of G. Note that by construction, F{As) = As holds if and only if G{Q*g + 
As) = GiQs) = 0. 

We now claim that F{M{r)) C B(r). Since F is continuous and B(r) is convex and compact, this inclusion 
implies, by Brouwer's fixed point theorem [20], that there exists some fixed point A5 G B(r). By uniqueness of the 
zero gradient condition (and hence fixed points of i^), we can thereby conclude that IIO5 — 0J||oo < r. 

Let A e MP^P denote the zero-padded matrix, equal to A5 on S* and zero on S'^. By definition, we have 

Gie*s + As) = -{e* + A)s^ + ^s + KZs 

= [ - (e* + A)s' + e*s] + [% - ms'] + xJs 

= [ - (e* + A)g' + e*s] +Ws + XnZs, (77) 

where we have used the definition W ~ — T,* . 
For any A5 G B(r), we have 



|||e*-^A||u < |||e*-^||UII|A||u 

< As- rf|lA|U, (78) 

where || A||oo denotes the elementwise £cx3-norm (as opposed to the £00 -operator norm |||A|||oo), and the inequality 
follows since A has at most d non-zero entries per row/column. 

By the definition (61) of the radius r, and the assumed upper bound (59), we have ||A||oo < r < sITTd' 
that the results of Lemma 5 apply. By using the definition (52) of the remainder, taking the vectorized form of the 
expansion (57), and restricting to entries in S, we obtain the expansion 

Yec{{e* + Ay' -e*-^)g+ T*ss As = vcc((e*-'A)Ve*-')^. (79) 
Using this expansion (79) combined with the expression (77) for G, we have 



F(As) = -{r*ssy G{e*s + As) + A 



s 



ns) 'vcc{[(e*+A) '-e*~']^-Ws-XnZs}+As 
= {T*ssy'yec[{Q*~'AfjQ*-']^ - (L^ -'(W^s + A„l. 



Ti T2 

The second term is easy to deal with: using the definition A'r* = |||(rg_5)^^|||oo, we have ||T2j|tx) < /\r»(||VF|| 
A„) = r/2. It now remains to show that UTiHoo < r/2. We have 

llTilU < Kr,\\vec[{Q*~'A)'jQ*-']J^ 
< Xr-P(A)||co, 

where we used the expanded form (57) of the remainder. Applying the bound (58) from Lemma 5, we obtain 



Since r < ^^■^ — ^ by assumption (59), we conclude that 

1 



3KI,.Kr^d 

thereby establishing the claim. 
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D Proof of Lemma 1 



For each pair and i' > 0, define tlie event 

1 " 

k=l 

As the sub-Gaussian assumption is imposed on the variables {X^''^} directly, as in Lemma A. 3 of Bickel and Levina 
[2], our proof proceeds by first decoupling the products x'>''^X^''\ For each pair (i, j), we define p*j — I]*^/^S*-S]*^-, 

and the rescaled random variables x'>''^ := X^j^ / -^S*-. Noting that the strict positive definiteness of S* implies that 
\p*j I < 1, we can also define the auxiliary random variables 



:=xi'=)+X?) and v!^^ := X^"^ - X?'^ 



(80) 



With this notation, we then claim: 

Lemma 9. Suppose that each X^'^'^ is sub-Gaussian with parameter a. Then for each node pair the following 
properties hold: 

(k) (k) 

(a) For all k = 1, . . . , n, the random variables U^j and V^j are sub-Gaussian with parameters 2a. 

(b) For all v > 0, the probability f[Kij{v)\ is upper bounded by 



p[iE(^?')'-2(i-p:,)i> 



2ni> 



k=l 



Proof, (a) For any r e M, we have 

E[exp(r[/J'^)] = E[exp(rXf^) exp (rxf ^) 



]+p[|^(yW)2_2(l_p*^.)l> 



2711^ 



k=l 



< E 



exp(2rXp^) 



1/2 r 



exp(2rXf^] 



1/2 



— (k) — (k) 

where we have used the Cauchy-Schwarz inequality. Since the variables X- and Xj ' are sub-Gaussian with param- 
eter cr, we have 



E 



exp (2rXf^) 



1/2 



cxp {2rX 



{k)^ 



1/2 



< cxp(a2_) cxp{cT^'—), 



(k) 

so that t/jj is sub-Gaussian with parameter 2a as claimed, 
(b) By straightforward algebra, we have the decomposition 



n 1 1 ^ 



fc=l i=l 



By union bound, we obtain that P[A.y (z^)] is upper bounded by 



n 

[l5:([/Wf-2(i+p:,.)|> 



Aril' 



fc=i 



n 

[lE(^f)'-2(i-p:,)|> 



which completes the proof of Lemma 9(b). 



It remains to control the terms J2^=ii^ij'^ )^ '^^^ Sfc=i C^if' ^'^ ^ly exploiting tail bounds [6] for sub- 

exponential random variables. A zero-mean random variable Z is said to be sub-exponential if there exists a constant 
7 G (0, oo) and G (0, oo] such that 



k=l 



.(fc)A2 



(81) 



□ 



E[exp(tZ)] < exp{-f^ t'^ /2) foralUG(-0,( 



(82) 
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Note that for (j) < +00, this requirement is a weakening of sub-Gaussianity, since the inequality is only required to 
hold on the interval (—0, +0). 

Now consider the variates Zk-^ij := [U^j^'Y — 2 (1 + p*^). Note that they are zero-mean; we also claim they are 
sub-exponential. 

Lemma 10. For a/Z fc £ {1, . . . , n} and node-pairs {i, j) V x V, the variables 

are sub-exponential with parameter = 16(1 + 4(7^) in the interval {—ipu, (j>u), with ~ 1/(16(1 + 4(7^)). 

We can exploit this lemma to apply tail bounds for sums of i.i.d. sub-exponential variates (Thm. 5.1, [6]). Doing 
so yields that for t < jI(I>u, we have P lELi(t^if^)^ - 2(1 + p*^)\ > nt] < 2exp{ - ni^/{2-fl)}. Setting 
t = 2u/ maxj S*-, and noting that (2?!;^/ > {2nh'/ max, S*-), we obtain 



\ym'r-2{i+p*,)\> 



k=l 



2nv 



< 2 cxp 



< 2exp 



2n,i/2 



max,(E*J2^2 



max,(Sl)2i28(l + 4(T2)2 



for all 1/ < 8(maxi S* ) (1 + 40-2). A similar argument yields the same tail bound for the deviation involving 
Consequently, using Lemma 9(b), we conclude that 

,2 



(fe) 
y ■ 



'[KjH] < 4cxp{ 



nv 



max,(S*)2l28(l + 4cr2)2 
vaUd for ly < 8(maxi S*j) (1 + 4cr2)^ as required. It only remains to prove Lemma 10 
Proof of Lemma 10. If we can obtain a bound _B > such that 



sup 

m>2 



1/rn 



it then follows (Thm. 3.2, [6]) that Zk-^ij is sub-exponential with parameter 2B in the interval We obtain 

such a bound B as follows. Using the inequality (a + 6)'" < 2™(a'" + 6™), valid for any real numbers a, b, we have 



< 2" 



(E(|C/^)|2™) + (2(l + p^)r)- 



(83) 



(k) 

Recalling that U^^ is sub-Gaussian with parameter 2(7, from Lemma 1 .4 of Buldygin and Kozachenko [6] regarding 
the moments of sub-Gaussian variates, we have E[|;7^'''|2™] < 2(2TO/e)"(2cr)2". Thus, noting the inequality m! > 
(m/e)™, it follows that E[\U^f'> |2™]/77i! < 2^"'+^a'^"'. It then follows from equation (83) that 



(m!)!/" 

where we have used the inequality {x + yY^™ < 2^/™(.t^/™ + t/^^'"), valid for any integer to e N and real numbers 
X, 1/ > 0. Since the bound is a decreasing function of m, it follows that 

1/™ / 4(1 + 



1/™ / 4('i I 

~ \ (m!)i/™ 

< 2i/™ ( 21/™16(t2 ' + 



sup 

m>2 



k;ij I 



< 2i/M2i/2l6a 



(2)1/2 



< 32cr2 + 8 = 8(1 + 4cr2), 



where we have used the fact that \p*j \ < 1. The claim of the lemma thus follows. 



□ 
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E Proof of Lemma 2 



(k) (k) (k) 

Define the random variables W^-^ = X- Xj — Zl*^ , and note that they have mean zero. By applying the Chebyshev 
inequality, we obtain 

n n 



r(fc) 

k=l fe=l 



< ^2LJL^ (84) 

Letting A — {(ai, . . . , a„) | G {0, . . . , 2m}, X)"^! ~ 2m}, by the multinomial theorem, we have 



fe=i neA ^ 1' ■ ■ ■ ' "/ /c=i 



^ai, . . . , a„ 



where the final equality uses linearity of expectation, and the independence of the variables {W^^^}^^l. 

Since the variables wj^^ are zero-mean, the product HLi ^[(Wi^/^)'"'] vanishes for any multi-index a ^ A such 



that flfe = 1 for at least one k. Accordingly, defining the subset 



A-i := {(ai, . . . ,a„) | G {0,2, . . .,2m}, ^ = 2?7i}, 



we have 



fc— 1 aGA-i A:— 1 

< ( 5: L J) max n (85) 



Ti T2 

The quantity Ti is equal to the number of ways to put 2m balls in n bins such that if a bin contains a ball, it should 
have at least two balls. Note that this implies there can then be at most m bins containing a ball. Consequently, the 
term Ti is bounded above by the product of the number of ways in which we can choose m out of n bins, and the 
number of ways in which we can put 2m balls into m bins — viz. 

Ti < rM m^™ < n'^m^'". 
\m J 

Turning now to the second term T2, note the following inequality: for any numbers {vi, . . . ,ve) G and non- 
negative integers (ai, . . . , ae), we have 

Ylvl" < ( max < Xl^^^' where a+ := flfe- 



fc=i 

k=l ' k=l 



k=l 
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Using this inequality, for any a G A-i, we have 

n 

k=l k=l 



{fc|afc#0} 

< m max E[(W,^V"]: 

where the last inequality follows since any multi-index a G A-i has at most m non-zero entries. Thus, we have shown 
thatTa < 77imaxfce{i,...,„}E[(VP^/f 

Substituting our bounds on Ti and T2 into equation (85), we obtain 

n 



mY.iK) 1 ^ n^m^'^+i max E[(W^,',"0""] ■ (86) 

fce{i,...,«} 



It thus remains to bound the moments of Wh ■ We have 



where we have used the inequality (a + 6)^™ < 2^'"(a^'" + 6^'"), valid for all real numbers a and b. An application 
of the Cauchy-Schwarz inequality yields 



E[(l^f)2™] < 22™{JE[(xfV"]E[(xf))'"'] + [sy2™} 



< 22™(i^„[S]^s^]™ + [S]:^f™) 



where we have used the assumed moment bound KUX^'^^/ VE^)^™] < Km- Equation (86) thus reduces to 



n 



k=l 

Substituting back into equation (84) yields 



1 " 



> iy\ < — — — 



fe=l 

Noting that E*^ , I]*^ , and E*^ are all bounded above by maxi E*;, we obtain 

{m2™+i22™(max,; E* )2™ (ii:„, + 1)} 



1 



I n 



> 1/ < 



as claimed. 
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